/* -----------------------------------------------------------------------
   See COPYRIGHT.TXT and LICENSE.TXT for copyright and license information
   ----------------------------------------------------------------------- */
/* -------------------------------------------------------------------------
    REFS:
    http://en.wikipedia.org/wiki/B-spline
    http://www.cs.mtu.edu/~shene/COURSES/cs3621/NOTES/surface/bspline-construct.html
    http://graphics.idav.ucdavis.edu/education/CAGDNotes/Quadratic-B-Spline-Surface-Refinement/Quadratic-B-Spline-Surface-Refinement.html
    ----------------------------------------------------------------------- */
#include "plm_config.h"
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "bspline.h"
#include "bspline_optimize.h"
#include "bspline_opts.h"
#if (CUDA_FOUND)
#include "bspline_cuda.h"
#include "cuda_util.h"
#include "delayload.h"
#endif
#include "mha_io.h"
#include "vf.h"
#include "fdk_opts.h"
#include "fdk_opts_a2p.h"
#include "drr_opts.h"
#include "proj_image_dir.h"
#include "bspline3DT2D_optimize.h"
#include "print_and_exit.h"
#include "proj_image.h"
#include "proj_image_set.h"
#include "fdk_cuda.h"
#include "drr_cuda.h"
#include "drr_util.h"
#include "fdk_util.h"
#include "hnd_io.h"
#include "read_projangle.h"
#include "intensity2att.h"
#include "marginpadding.h"
#include "ramp_filter.h"
#include "zhang_cuda.h"
#include "pvalue_optimize.h"
#include "proj_image.h"
#include <iostream>
#ifndef _WIN32
#include <dlfcn.h>
#endif

#define HALFCENTER 64
#define MARGIN 5

int
main (int argc, char* argv[])
{
	Fdk_options cb_options;	
	Drr_options drr_options;
	Proj_image_dir *proj_dir;
	void * fdk_dev_state;
	void * drr_dev_state;

    Bspline_options options;
    Bspline_parms *parms = &options.parms;
    Bspline_xform *bxf;
    Volume *moving, *fixed, *moving_grad;
    Volume *vector_field = 0;
    Volume *moving_warped = 0;
    int roi_offset[3];

	printf("parsing args\n");

	int last_arg=fdk_parse_args_find_last(&cb_options, argc, argv);
    bspline_opts_parse_args (&options, argc-last_arg, argv+last_arg);

    

	//copy bspline threading to fdk threading
	if (options.parms.threading==BTHR_CUDA)
	cb_options.threading=THREADING_CUDA;
	else
	cb_options.threading=THREADING_CPU_OPENMP;
	

#if (CUDA_FOUND)
    if (parms->threading == BTHR_CUDA) {
        if (!delayload_cuda()) { exit(0); }
        LOAD_LIBRARY (libplmcuda);
        LOAD_SYMBOL (CUDA_selectgpu, libplmcuda);
        CUDA_selectgpu (parms->gpuid);
        UNLOAD_LIBRARY (libplmcuda);
    }
#endif

    fixed = read_mha (options.fixed_fn);
    if (!fixed) exit (-1);
    moving = read_mha (options.moving_fn);
    if (!moving) exit (-1);

    volume_convert_to_float (moving);
    volume_convert_to_float (fixed);

	Volume * vol = my_create_volume (&cb_options);
	Volume * vol2 = my_create_volume (&cb_options);

    printf ("Making gradient\n");
    moving_grad = volume_make_gradient (moving);

#if defined (commentout)
    /* Load and adjust landmarks */
    if (options.fixed_landmarks && options.moving_landmarks) {
	parms->landmarks = bspline_landmarks_load (
	    options.fixed_landmarks, options.moving_landmarks);
	bspline_landmarks_adjust (parms->landmarks, fixed, moving);
    }
#endif

    /* Debug */
    //write_mha ("moving_grad.mha", moving_grad);

    /* Allocate memory and build lookup tables */
	printf("flag\n\n\n");
    printf ("Allocating lookup tables\n");
    memset (roi_offset, 0, 3*sizeof(int));
    if (options.input_xf_fn) {
	bxf = bspline_xform_load (options.input_xf_fn);
	if (!bxf) {
	    fprintf (stderr, "Failed to load %s\n", options.input_xf_fn);
	    exit (-1);
	}
    } else {
	bxf = (Bspline_xform*) malloc (sizeof (Bspline_xform));
	bspline_xform_initialize (
	    bxf,
	    fixed->offset,
	    fixed->spacing,
	    fixed->dim,
	    roi_offset,
	    fixed->dim,
	    options.vox_per_rgn
	);
    }



	/* Look for input files */
	proj_dir = proj_image_dir_create (cb_options.input_dir);
	if (!proj_dir) {
		print_and_exit ("Error: couldn't find input files in directory %s\n",
			cb_options.input_dir);
	}

	/* Choose subset of input files if requested */
	if (cb_options.image_range_requested) {
		proj_image_dir_select (proj_dir, cb_options.first_img, 
			cb_options.skip_img, cb_options.last_img);
	}


	float angles[1000];
	double sad=1000.0;
	double sid=1500.0;

	char projfilename[128];
	strcpy(projfilename,cb_options.input_dir);
	strcat(projfilename,"\\ProjAngles.txt");
	int nprj=read_ProjAngle(projfilename,angles);
	// This is just to retrieve the 2D image dimensions
	Proj_image *cbitmp = proj_image_dir_load_image (proj_dir, 0);
	cbitmp->pmat  = proj_matrix_create ();	




	CB_Set PSd;
	CB_Set PSd2;


	LOAD_LIBRARY (libplmcuda);
	LOAD_SYMBOL (fdk_cuda_state_create, libplmcuda);
	LOAD_SYMBOL (fdk_cuda_queue_image, libplmcuda);
	LOAD_SYMBOL (fdk_cuda_fetch_volume, libplmcuda);
	LOAD_SYMBOL (fdk_cuda_backproject, libplmcuda);
	LOAD_SYMBOL (fdk_cuda_state_destroy, libplmcuda);

	fdk_dev_state = fdk_cuda_state_create (
		vol, 
		cbitmp->dim[0]*cbitmp->dim[1], 
		1, 
		&cb_options
		);

	char img_fn[256];

	drr_opts_init (&drr_options);
	set_image_parms (&drr_options);

	PSd.nproj=0;
	int angle_i=cb_options.first_img;
	for (int i=0; i<proj_dir->num_proj_images; i++)
	{
		PSd.cbi[i]=proj_image_create();
		raw_load(PSd.cbi[i],proj_dir->proj_image_list[i]); 

		marginpadding(PSd.cbi[i]);
		intensity2att(PSd.cbi[i]);


		//sprintf (img_fn, "%s%04d.raw", drr_options.output_prefix,680-i );
		//raw_save (PSd.cbi[i], img_fn);
		PSd.nproj++;

		//Half-fan conversion
		if(!cb_options.full_fan)
		{
			
			float buffer[512*384];
			memcpy(buffer,PSd.cbi[i]->img,sizeof(float)*512*384);

			PSd.cbi[i]->dim[0]=512-HALFCENTER;
	
			for(int r=0;r<PSd.cbi[i]->dim[1];r++)
				for(int c=0; c<PSd.cbi[i]->dim[0]; c++)
					(PSd.cbi[i]->img)[r*PSd.cbi[i]->dim[0]+c]=(buffer)[r*cbitmp->dim[0]+c+HALFCENTER];
	
		}

		//// applied ramp-filter
		RampFilter_filter_extend_DC(PSd.cbi[i]->img,PSd.cbi[i]->img,PSd.cbi[i]->dim[0],PSd.cbi[i]->dim[1],PSd.nproj);

		//char * mat_fn=0;
		//proj_image_save (PSd.cbi[0], "C:\\AAAFiles\\CBCTData\\Proj.raw", mat_fn);
		
		PSd.cbi[i]->pmat = proj_matrix_create ();	
		set_proj_matrix (PSd.cbi[i], (double)angles[angle_i], sad, sid, cb_options.full_fan);

		angle_i+=cb_options.skip_img;


	}

	
	CB_Set PSc;
	CB_Set_clone(&PSc, &PSd);


	//CUDA_Bkpj(vol,&PSd,fdk_dev_state);

	//write_mha("C:\\AAAFiles\\CBCTData\\CBCT2.mh5",vol);



	//memset(vol->img,0,sizeof(float)*vol->npix);
	
	//vol=read_mha("C:\\AAAFiles\\CBCTData\\CBCT1.mh5");	
	//vol2=read_mha("C:\\AAAFiles\\CBCTData\\CBCT2.mh5");	

	//for (int i=0; i<vol->npix; i++){
	//	((float*)vol->img)[i]=0.0f;
	//	((float*)vol2->img)[i]=0.0f;
	//}

	//((float*)vol->img)[128*128*50+128*64+64]=1.0f;
	//((float*)vol2->img)[128*128*50+128*64+64]=1.0f;


	////half fan
	drr_options.image_center[0]=0;
	drr_options.image_center[1]=0;


	//memset(vol->img, 0, vol->npix*sizeof(float));
	//for (int z=0; z<vol->dim[2]; z++)
	//	for (int y=0; y<vol->dim[1]; y++)
	//		for (int x=0; x<vol->dim[0]; x++){
	//			if ((x-vol->dim[0]/2)*(x-vol->dim[0]/2)+(y-vol->dim[0]/2)*(y-vol->dim[0]/2)<vol->dim[0]/4*vol->dim[0]/4)
	//				((float *)vol->img)[z*vol->dim[0]*vol->dim[1]+y*vol->dim[0]+x]=1;
	//		}

	//write_mha("C:\\AAAFiles\\CBCTData\\CBCT0.mh5",vol);

	//commmented out temporarily

	//bool copyvol2dev=true;

	//drr_dev_state =drr_cuda_state_create_cu_from_fdk(fdk_dev_state,vol,copyvol2dev);

	//CUDA_Pj2(&PSc,vol, drr_dev_state);
	//CUDA_Pj2(&PSd,vol2, drr_dev_state);


	//double PSc_by_PSd=CB_Set_Prod(&PSc,&PSd);

	//CUDA_Bkpj(vol2,&PSd,fdk_dev_state);


	CUDA_Bkpj(vol2,&PSd,fdk_dev_state);

	write_mha("C:\\AAAFiles\\CBCTData\\CBCT2BP.mh5",vol2);

	double V1_by_V2=Volume_Prod(vol,vol2);

	//if (!delayload_cuda ()) { exit (0); }
	//fdk_cuda_fetch_volume (dev_state, vol->img, vol->npix * vol->pix_size);


	/* Free memory on device */
	//if (!delayload_cuda ()) { exit (0); }
	//fdk_cuda_state_destroy (fdk_dev_state);

	//proj_image_destroy(cbitmp);

	//write_mha("C:\\AAAFiles\\CBCTData\\CBCT0.mh5",vol);



    /* Run the optimization */

	printf ("Running optimization.\n");
//    bspline3DT2D_optimize (bxf, 0, parms, fixed, moving, moving_grad, &cb_options);

//	CUDA_Pj(&PSd,vol, drr_dev_state, &drr_options);

	pvalue_optimize (bxf, 0, parms, vol, &PSd, moving, moving_grad, fdk_dev_state, drr_dev_state, &cb_options, &drr_options);


    printf ("Done running optimization.\n");

    /* Save output transform */
    if (options.output_xf_fn) {
	bspline_xform_save (bxf, options.output_xf_fn);
    }

    /* Create vector field from bspline coefficients and save */
    if (options.output_vf_fn 
	|| options.output_warped_fn 
#if defined (commentout)
	|| (options.warped_landmarks && options.fixed_landmarks 
	    && options.moving_landmarks)
#endif
    )
    {
	printf ("Creating vector field.\n");
	vector_field = new Volume (fixed->dim, fixed->offset, 
	    fixed->spacing, fixed->direction_cosines, 
	    PT_VF_FLOAT_INTERLEAVED, 3);
	if (parms->threading == BTHR_CUDA) {
#if (CUDA_FOUND)
	    LOAD_LIBRARY (libplmcuda);
	    LOAD_SYMBOL (CUDA_bspline_interpolate_vf, libplmcuda);
	    CUDA_bspline_interpolate_vf (vector_field, bxf);
	    UNLOAD_LIBRARY (libplmcuda);
#else
	    bspline_interpolate_vf (vector_field, bxf);
#endif
	} else {
	    bspline_interpolate_vf (vector_field, bxf);
	}
    }

    /* Create warped output image and save */
    if (options.output_warped_fn) {
	printf ("Warping image.\n");
	moving_warped = vf_warp (0, moving, vector_field);
	if (moving_warped) {
	    printf ("Writing warped image.\n");
	    write_mha (options.output_warped_fn, moving_warped);
	} else {
	    printf ("Sorry, couldn't create warped image.\n");
	}
    }

    /* Write the vector field */
    if (options.output_vf_fn) {
	printf ("Writing vector field.\n");
	write_mha (options.output_vf_fn, vector_field);
    }

    /* Free memory */
    printf ("Done warping images.\n");
    bspline_parms_free (parms);
    bspline_xform_free (bxf);
    free (bxf);
    delete fixed;
    delete moving;
    delete moving_grad;
    delete moving_warped;
    delete vector_field;

    printf ("Done freeing memory\n");

    return 0;
}
